In [1]:
import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import sys
import cPickle
import scipy.io as sio
import matplotlib.units
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
from common_functions import *
import half_nanoplate_functions as hnf
import matplotlib_helper_functions as mplhf
%matplotlib inline
#%config InlineBackend.figure_format = 'retina'
#import seaborn as sns
C:\Users\Scherer Lab E\Anaconda2\envs\170112\lib\site-packages\skimage\viewer\utils\core.py:10: UserWarning: Recommended matplotlib backend is `Agg` for full skimage.viewer functionality.
  warn("Recommended matplotlib backend is `Agg` for full "
In [2]:
os.chdir("C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Data")
In [3]:
import cPickle


f = open("force_distributions_expt.pkl", "r")
force_dist = cPickle.load(f)
f.close()

f = open("log_pd_experiment.pkl", "r")
log_pd_ls, edge_dist_ls, theta_ls, edge_theta_ls, avg_r_ls = cPickle.load(f)
f.close()

store = pd.HDFStore("K:\Pat's_Projects\ParticleTrajectoryData\half_nanoplate_dynamics_processed.h5", mode='r')
keys = store.index['key']
np_pos = [store.get(v).drop_duplicates(['frame', 'track id']).copy() for v in keys[:-1]]
store.close()

pmf_345 = mplhf.load_matlab_to_dict("pmf_barrier_kernelfit_L345.mat")

umb_samp_pmf = pd.read_csv("emus_pmfs.txt", sep=' ')

Derivatives of Umbrella Smapled pmfs

In [146]:
def calculated_derivative(x, y, point_shift=7):
    """Calculates the derivative of two Series where one is x and the
    other is f(x).
    
    :param x: A pandas Series object of a independent variable
    :param y: A pandas Series object of a dependent variable such that
    y = f(x).
    :param point_shift: The number that points will be shifted for the 
    derivative. At point_shift=1 the instantanous change in x and y are 
    calculated. At point_shift=3 the change of x and y are caclulated 3
    points in the future."""
    
    if (type(x).__module__, type(y).__module__) == (np.__name__, np.__name__):
        change_x = x[point_shift:] - x[:-point_shift]
        change_y = y[point_shift:] - y[:-point_shift]
        x_points = x[:-point_shift] + (change_x[point_shift+1]/2.0)
    else:
        change_x = x.shift(point_shift) - x
        change_y = y.shift(point_shift) - y
        x_points = x + (change_x[point_shift+1]/2.0)
    derivative = change_y/change_x
    
    
    return x_points, derivative

Derivative of Umbrella Sampled pmfs

In [152]:
results = []
for l in range(6):
    results_dict = {}
    
    # Clean up the US data
    x_us = umb_samp_pmf.y*10**3 - 937.5
    pmf_us = umb_samp_pmf[str(l)] / (4.11*10**-21)
    
    x_1st_der, der_1st = calculated_derivative(x_us, pmf_us, point_shift=7)
    x_2nd_der, der_2nd = calculated_derivative(x_1st_der, der_1st, point_shift=7)
    
    x_spline = np.linspace(-600, 600, 2000)
    spline_prams = scipy.interpolate.splrep(x_us, pmf_us, s=25)
    spline_y = scipy.interpolate.splev(x_spline, spline_prams)
    spline_der_1st = scipy.interpolate.splev(x_spline, spline_prams, der=1)
    spline_der_2nd = scipy.interpolate.splev(x_spline, spline_prams, der=2)
    
    results_dict['x'] = x_us
    results_dict['pmf'] = pmf_us
    results_dict['x_1st_der'] = x_1st_der
    results_dict['x_2nd_der'] = x_2nd_der
    results_dict['1st_der'] = der_1st
    results_dict['2nd_der'] = der_2nd
    results_dict['x_spline'] = x_spline
    results_dict['spline'] = spline_y
    results_dict['1st_der_spline'] = spline_der_1st
    results_dict['2nd_der_spline'] = spline_der_2nd
    
    results.append(results_dict)
    
    plt.figure(figsize=(4, 6))
    plt.subplot(311)
    plt.plot(x_us, pmf_us)
    plt.plot(x_spline, spline_y)
    plt.title('L='+str(l))
    plt.xlabel('Distance from Edge (nm)')
    plt.ylabel(r'pmf (kbt)')
    
    plt.subplot(312)
    plt.plot(x_1st_der, der_1st, 'o')
    plt.plot(x_spline, spline_der_1st)
    plt.plot((-600, 600), (0,0), '--k')
    plt.xlabel('Distance from Edge (nm)')
    plt.ylabel(r'$\frac{d}{dx}\mathregular{pmf_{US}}$')
    
    plt.subplot(313)
    plt.plot(x_2nd_der, der_2nd, 'o')
    plt.plot(x_spline, spline_der_2nd)
    plt.plot((-600, 600), (0,0), '--k')
    plt.xlabel('Distance from Edge (nm)')
    plt.ylabel(r'$\frac{d^2}{dx^2}\mathregular{pmf_{US}}$')
    
    plt.show()
In [130]:
def nearest_index_in_array(val, array):
    return np.argmin(np.abs(array - val))

Calculations from 1st Derivative

Calculate where the first derivative crosses 0 for L=0-3 to find the minima and maxima that make up the barrier. From this can calculate the $\Delta U$ and $x_b$ from each pmf. The derivative is calculated from the data directly. Zero crossing point in x is found using linear interpolation between the two points that are nearest to crossing 0. The x values are used in a spline fit of the original data to calculate $\Delta U$.

In [131]:
x_bs = []
del_gs = []

for num, l_dict in enumerate(results):
    deriv = l_dict['1st_der']
    x_deriv = l_dict['x_1st_der']
    
    
    zero_crossing = deriv[(deriv.shift(1) * deriv) < 0]
    if len(zero_crossing) <= 1:
        continue
    index = zero_crossing[:2].index
    x_crossing = []
    
    # Find where it crosses 0 using linear interpolation
    for idx in index:
        f = scipy.interpolate.interp1d(deriv.loc[idx-1:idx].values, x_deriv.loc[idx-1:idx].values, kind='linear')
        x_cross = f(0)
        x_crossing.append(x_cross)
    
    original_indx = [nearest_index_in_array(v, l_dict['x_spline']) for v in x_crossing]
    
    x_b = x_crossing[-1] - x_crossing[0]
    org_pmf = l_dict['spline']
    del_g = org_pmf[original_indx[1]] - org_pmf[original_indx[0]]
    
    
    print "L = "+str(num)
    print "Minimum at "+str(x_crossing[0])+" nm and "+str(org_pmf[original_indx[0]])+" kbT"
    print "Maximum at "+str(x_crossing[1])+" nm and "+str(org_pmf[original_indx[1]])+" kbT"
    
    
    print "x_b = "+str(x_b)
    print "del_g = "+str(del_g)+"\n"
    
    x_bs.append(x_b)
    del_gs.append(del_g)
L = 0
Minimum at -498.04009457 nm and 74.0989647892 kbT
Maximum at -83.3823520668 nm and 100.737354051 kbT
x_b = 414.657742503
del_g = 26.6383892615

L = 1
Minimum at -429.597211902 nm and 110.052297559 kbT
Maximum at -122.699006059 nm and 122.012101703 kbT
x_b = 306.898205843
del_g = 11.9598041442

L = 2
Minimum at -272.460363134 nm and 131.014814885 kbT
Maximum at -162.812382889 nm and 132.742591664 kbT
x_b = 109.647980244
del_g = 1.72777677962

L = 3
Minimum at -224.170684394 nm and 137.636186353 kbT
Maximum at -217.339281736 nm and 137.587642191 kbT
x_b = 6.83140265786
del_g = -0.048544162934

Calculations from 2nd Derivative

The minima and the maxima in the 2nd derivative are found as they should correspond to the minima and maxima in the original function. A maxima in the 2nd derivative corresponds to a minimum in the original function and a minimum in the 2nd derivative corresponds to a maximum in the original function. Some of these second derivative curves have two maxima near where a minimum is expected in the orginal function. The one chosen is the one with the highest x value. No interpolation was done in this analysis as the max/min was just found on an array (within some bounds) to determine the maxima and minima of the 2nd derivative. From the x values of the maxima and minima, $x_b$ and $\Delta U$ were calculated.

In [138]:
x_bs_2nd = []
del_gs_2nd = []

for num, l_dict in enumerate(results):
    x = l_dict['x_2nd_der']
    y = l_dict['2nd_der']
    maxima = np.argmax(y[np.logical_and(x < -200, x > -350)])
    if num <= 1:
        maxima = np.argmax(y[np.logical_and(x < -200, x > -500)])
    
    
    minima = np.argmin(y[np.logical_and(x < 0, x > -100)])
    
    indx_min = nearest_index_in_array(x[minima], l_dict['x'])
    indx_max = nearest_index_in_array(x[maxima], l_dict['x'])
    
    pmf = l_dict['pmf']
    
    x_b = x[minima] - x[maxima]
    del_g = pmf[indx_min] - pmf[indx_max]
    
    print "\nL="+str(num)
    print 'x_b = '+str(x[minima] - x[maxima])
    print 'del_g = '+str(pmf[indx_min] - pmf[indx_max])
    
    x_bs_2nd.append(x_b)
    del_gs_2nd.append(del_g)
L=0
x_b = 377.0
del_g = 24.1639392708

L=1
x_b = 377.0
del_g = 10.8556157079

L=2
x_b = 234.0
del_g = -1.41109855649

L=3
x_b = 208.0
del_g = -3.25642848101

L=4
x_b = 221.0
del_g = -15.8335982968

L=5
x_b = 169.0
del_g = -14.3206943865

Plot $\Delta G^\ddagger$ and $x^\ddagger$ vs $l$ from 1st Derivative

Taken from the first derivative analysis.

In [143]:
fig = plt.figure(figsize=(3.35, 2))
Ls = range(len(x_bs))
ax1 = plt.subplot()
plt.plot(Ls, del_gs)
plt.xlabel(r'$l$', usetex=True)
plt.ylabel(r'$\Delta G^\ddagger$ or $\Delta U$ ($k_B T$)', usetex=True)

ax1_2 = ax1.twinx()
ax1_2.plot(Ls, x_bs, 'r')
plt.ylabel(r'$x^\ddagger$ (nm)', usetex=True, fontdict={'color':'red'})

plt.show()

params = {}
params['ls'] = Ls
params['x_bs'] = x_bs
params['del_gs'] = del_gs
In [144]:
'''Save the componenets of the del_g and x_0 graph'''
import cPickle

pickle_file = open("C:\Users\Scherer Lab E\Box Sync\Half-Nanoplate\Figures and Data\Data\\graph_del_g_x_o.pkl", 'w')
cPickle.dump(params, pickle_file)
pickle_file.close()

Plot $\Delta G^\ddagger$ and $x^\ddagger$ vs $l$ from 2nd Derivative

Taken from the second derivative analysis.

In [141]:
fig = plt.figure(figsize=(3.35, 2))
Ls = range(len(x_bs_2nd))
ax1 = plt.subplot()
plt.plot(Ls, del_gs_2nd)
plt.xlabel(r'$l$', usetex=True)
plt.ylabel(r'$\Delta G^\ddagger$ or $\Delta U$ ($k_B T$)', usetex=True)

ax1_2 = ax1.twinx()
ax1_2.plot(Ls, x_bs_2nd, 'r')
plt.ylabel(r'$x^\ddagger$ (nm)', usetex=True, fontdict={'color':'red'})

plt.show()